home *** CD-ROM | disk | FTP | other *** search
/ QRZ! Ham Radio 8 / QRZ Ham Radio Callsign Database - Volume 8.iso / pc / files / dsp / 56000tar.z / 56000tar / 56000 / speech / durbin1.asm < prev    next >
Assembly Source File  |  1991-11-26  |  6KB  |  156 lines
  1. ;
  2. ; This program originally available on the Motorola DSP bulletin board.
  3. ; It is provided under a DISCLAMER OF WARRANTY available from
  4. ; Motorola DSP Operation, 6501 Wm. Cannon Drive W., Austin, Tx., 78735.
  5. ; 
  6. ; Durbin Solution for PARCOR (LPC) Coefficients
  7. ; 
  8. ; Last Update 3 June 87   Version 1.2
  9. ;
  10. ;       Implements the Durbin LPC algorithm.
  11. ;
  12.         opt     cex
  13.         page    132,66,3,3
  14. ;
  15. ;       i/o file definitions
  16. ;
  17. datain   equ    $ffff
  18. dataout  equ    $fffe
  19. nk       equ    10              ;10th order LPC analysis
  20.  
  21.         org     x:0             ;I/O values in X: memory
  22. r       ds      nk+1            ;autocorrelation coefficients
  23. k       ds      nk              ;reflection coefficients
  24.  
  25.         org     y:0             ;temporary values in Y: memory
  26. acoeffs dc      $7fffff         ;define a[0] = 1.0
  27.         ds      nk              ;LPC filter coefficients
  28. anew    dc      $7fffff         ;define anew[0] = 1.0
  29.         ds      nk              ;updated LPC filter coefficients
  30. alpha   ds      1               ;Durbin alpha
  31. error   ds      1               ;Durbin error
  32.  
  33.         org     p:$100
  34. ;
  35. ;       find out how many frames to process
  36. ;
  37.         movep   y:datain,n0
  38.         do      n0,endframe
  39.  
  40. ;
  41. ;       load the autocorrelations from an external file
  42. ;
  43.         move    #r,r0           ;load autocorrelation
  44.         do      #nk+1,loadr     ;get nk+1 autocorrelation values
  45.         movep   y:datain,a      ;get next autocorrelation value
  46.         move    a,x:(r0)+       ;and save it in the r array
  47. ;
  48. ;       The following registers correspond to the C variables:
  49. ;
  50. ;       r0 - r[i]               r4 - acoeffs[i]
  51. ;       r1 - k[i-1]             r5 - acoeffs[j-1]
  52. ;       r2 - r[i-j+1]           r6 - acoeffs[i-j+1]
  53. ;       r3 - anew[j-1]          r7 - loop counter
  54. ;
  55. ;       Approximate analysis times:
  56. ;
  57. ;            8th  order: 63.5  uS
  58. ;            10th order: 85.6  uS
  59. ;            16th order: 165.7 uS
  60. ;
  61. ;       Initialization
  62. ;
  63. loadr   move    #r,r0                   ;r(0) points to r[0]
  64.         move    #k,r1                   ;r(1) points to k[0]
  65. ;
  66. ;       Begin Durbin Algorithm
  67. ;
  68.         move    x:(r0)+,x0              ;get r[0]
  69.         move    x:(r0)+,a               ;get r[1], r(0) points to r[i]
  70.         abs     a          a,b          ;get abs(r[1]), copy r[1] to b
  71.         eor     x0,b  #acoeffs+1,r4     ;N = sign bit, r(4) points to a[i]
  72.         and     #$fe,ccr                ;make quotient positive
  73.         rep     #24                     ;set up for 24-bit quotient
  74.         div     x0,a                    ;get abs(r[1])/r[0]
  75.         jmi     L1                      ;check sign bit N
  76.         neg     a                       ;negate quotient if needed
  77. L1      move    a0,a                    ;move -r[1]/r[0] to A
  78.         clr     b   a,x:(r1)+  a,y0     ;k[0] = -r[1]/r[0], copy k[0] to y0
  79.         move    #$800000,b1             ;b = 1.0, r(1) points to k[i-1]
  80.         macr    -y0,y0,b  #2,r7         ;b=1.0-(k[0]*k[0]),loop counter=2
  81.         move    b,x1  a,y:(r4)+         ;x1 = 1.0-(k[0] * k[0]), a[1] = k[0]
  82.         mpyr    x1,x0,a   #2,n7         ;a1 = r[0] * (1.0 - (k[0] * k[0]))
  83.         move    a,y1                    ;alpha = r[0] * (1.0 - k[0] * k[0]))
  84.         move    #-2,n5                   ;initialize n(5) 
  85. ;
  86. ;       outer do loop  (note: alpha = y1)
  87. ;
  88.         do      #nk-1,L6                ;do outer do loop (nk-1) times
  89.         move    r0,r2                   ;r(2) points to r[i]
  90.         move    #acoeffs,r5             ;r(5) points to a[0]
  91.         move    (r0)+                   ;r(0) points to next r[i]
  92.         clr     a x:(r2)-,x0 y:(r5)+,y0 ;error = 0, preload 1st operand set
  93. ;
  94. ;       inner do loop #1  (note: r7 = i)
  95. ;
  96.         do      r7,L2                         ;do inner do loop #1 (i) times
  97.         mac     x0,y0,a x:(r2)-,x0 y:(r5)+,y0 ;error += a[j-1] * r[i-j+1]
  98. ;
  99. ;       back to outer do loop  (note: error = a)
  100. ;
  101. L2      abs     a          a,b          ;get abs(error), copy error to b
  102.         eor     y1,b       #2,n6        ;N = sign bit, y1 = alpha
  103.         and     #$fe,ccr                ;make quotient positive
  104.         rep     #24                     ;set up for 24-bit quotient
  105.         div     y1,a                    ;get abs(error)/alpha
  106.         jmi     L3                      ;check sign bit N
  107.         neg     a                       ;negate quotient if needed
  108. L3      move    a0,a                    ;move k[i-1] = -(error/alpha) to a
  109.         clr b   a,x0    a,y:(r4)+       ;put k[i-1] in x0, a[i] = k[i-1]
  110.         move    a,x:(r1)+  y:(r7)-,a    ;k[i-1] = -(error/alpha), dec r(7)
  111.         move    #$800000,b1             ;b = 1.0
  112.         macr    -x0,x0,b    r4,r6       ;b = 1.0 - (k[i-1] * k[i-1])
  113.         move    b,x1                    ;r(6) points to a[i+1]
  114.         mpyr    x1,y1,b   (r6)-n6       ;b = alpha * (1.0-(k[i-1]*k[I-1]))
  115.         move    b,y1                    ;save alpha=alpha*(1-k[i-1]*k[i-1])
  116.         move    #acoeffs+1,r5           ;r(5) points to a[j-1]
  117.         move    #anew+1,r3              ;r(3) points to anew[j-1]
  118.         move    y:(r6)-,y0              ;y0 = a[i-j+1], x0 = k[i-1]
  119.         move    y:(r5)+,a               ;a = a[j-1]
  120. ;
  121. ;       inner do loop #2  (note: r7 = (i-1))
  122. ;
  123.         do      r7,L4                   ;do inner do loop #2 (i-1) times
  124.         macr    x0,y0,a  y:(r6)-,y0     ;get anew[j-1], y0 = next a[i-j+1]
  125.         move    a,y:(r3)+               ;anew[j-1]=a[j-1]+k[i-1]*a[i-j+1]
  126.         move    y:(r5)+,a               ;a = a[j-1]
  127. ;
  128. ;       end of inner do loop #2
  129. ;
  130. ;       inner do loop #3  (note: r7 = (i-1))
  131. ;
  132. L4      move    x:(r3)-,x0 y:(r5)+n5,y0 ;dummy reads to dec r(3) and r(5)
  133.         do      r7,L5                   ;do inner do loop #3 (i-1) times
  134.         move    y:(r3)-,a               ;get anew[j-1]
  135.         move    a,y:(r5)-               ;a[j-1] = anew[j-1]
  136. ;
  137. ;       end of inner do loop #3
  138. ;
  139. ;       end of outer do loop
  140. ;
  141. ;       end of analysis
  142. ;
  143. ;       output reflection coefficients to an external file
  144. ;
  145. L5      move    (r7)+n7                 ;update loop counter
  146. L6      move    #k,r1                   ;point to reflection coeffs.
  147.         do      #nk,endsend             ;output nk k values
  148.         move    x:(r1)+,x0              ;get next k value
  149.         move    x0,y:dataout            ;output next K coefficient value
  150. endsend
  151.         nop
  152. endframe
  153.         nop
  154.         end
  155.  
  156.